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Abstract. In the probabilistic approach to quantum many-body systems, the 
ground-state energy is the solution of a nonlinear scalar equation written either as a 
cumulant expansion or as an expectation with respect to a probability distribution 
of the potential and hopping (amplitude and phase) values recorded during an 
infinitely lengthy evolution. We introduce a perturbative expansion of this probability 
distribution which conserves, at any order, a multinomial-like structure, typical of 
uncorrelated systems, but includes, order by order, the statistical correlations provided 
by the cumulant expansion. The proposed perturbative scheme is successfully tested 
in the case of spin Vi hard-core boson Hubbard models. 
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1. Introduction 

A multitude of evolution problems, including quantum many-body systems, can be cast 
in the form of a linear flow, namely a system of linear differential equations with respect 
to a parameter, the time, which can be real or imaginary. The solution of a linear flow, 
with a real or an imaginary time, admits an exact probabilistic representation, namely 
a Feynman-Kac-like formula, in terms of a proper collection of independent Poisson 
processes [1, 2, 3, 4, 5]. For a lattice system, the Poisson processes are associated with 
the links of the lattice and the probabilistic representation leads to an optimal algorithm 
[3, 4, 5] which coincides with the Green function quantum Monte Carlo method in the 
limit when the latter becomes exact [6]. The algorithm can be rigorously generalized [7] 
to include fluctuation control techniques, like reconfigurations and importance sampling, 
and allows the exact simulation of time-dependent correlation functions for system not 
affected by the so called sign problem [8]. 

In the limit of an infinitely-long imaginary time, the above exact probabilistic 
representation has been developed to yield semi-analytical results. In fact, for an 
arbitrary many-body system we are able to relate the energy of its ground state to 
the unique solution of a nonlinear scalar equation [9, 10, 11, 12]. This equation can be 
written in terms of a series involving the cumulants of integers, the multiplicities Ny, 
and N\, which count how many times the potential, hopping and phase variables take 
the values V, T and A during an infinitely long evolution of the system. The potential 
variables are, in some chosen base, the diagonal matrix elements of the Hamiltonian 
of the system, whereas the hopping and phase variables are related to the amplitude 
and phase of the off-diagonal matrix elements. Alternatively, the equation for the 
ground-state energy can be written in terms of an expectation involving the probability 
distribution of the multiplicities Ny, and N\. 

The two ways of writing the equation for Eq, cumulant expansion or probabilistic 
expectation, correspond to two different approaches to evaluate the ground state of a 
many-body system. In the former case, we can imagine to measure the exact cumulants 
of the system up to some finite (small) order, insert them in a corresponding truncated 
equation and solve it obtaining an approximation to E [11]. In the latter case, we have 
the possibility to make some guess on the probability distribution of the multiplicities 
Ny, Nt and N\ bypassing the microscopic connection between configurations of the 
system and values of the variables V, T and A. The simplest guess is to neglect any 
correlation among the multiplicities and assume that they are multinomially distributed. 
There is a class of systems for which, in the thermodynamic limit, a multinomial 
distribution exactly applies. The class includes the uniformly fully connected models, 
namely a collection of states all connected with equal hopping coefficients and in the 
presence of a potential operator with arbitrary levels and degeneracies, and the random 
potential systems, in which the hopping operator is generic and arbitrary potential levels 
are assigned randomly to the states with arbitrary probabilities. For this class of models 
we have found a zero-temperature universal thermodynamic limit displaying a quantum 
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phase transition [12]. 

Both the two approaches described above have limitations. The more severe 
drawback of the cumulant expansion is that truncating the series corresponds to 
introducing a rather artificial equivalent system with zero cumulants of large order. Since 
the determination of E is a problem which involves large fluctuations [9], we expect, 
and find indeed, some nonphysical behavior of the solutions corresponding to large 
interaction energies. On the other hand, we expect that the multinomial probability 
distribution used in the second approach may give quantitatively inaccurate results for 
most of the systems, i.e. when the correlations among the potential, hopping and phase 
multiplicities cannot be neglected. 

In the present paper we introduce a perturbative scheme merging the merits of 
the cumulant expansion with those of the expectation taken from an uncorrelated 
multinomial distribution. The idea is to develop a perturbative expansion of the 
probability distribution of the multiplicities N v , N T and N\ which conserves, at any 
order, a multinomial-like structure. Order by order, we add correlations among the 
multiplicities in such a way to modify the cumulants of the multinomial-like distribution 
and make them to coincide with those measured in the system up to the order considered. 
As a result, we gain better and better approximations to the real probability distribution. 
At the first order, we have a probability distribution which contain infinitely many 
cumulants and the first one is exact; at the second order, we have a probability 
distribution with infinitely many cumulants and the first two are exact; and so on. In 
this way, we expect to obtain, even at very small perturbative orders, accurate results for 
the ground-state energy of both weakly- and strongly-interacting many-body systems. 
We have checked our perturbative scheme in the case of spin V2 hard core boson Hubbard 
models in two dimensional square lattices. It is remarkable that, already at the second 
perturbative order, we find a ground-state energy which compares rather well with the 
exact value of E . 

The main advantage of our method lies in its semi-analytical character. Once all 
the cumulants up to some order k are measured, via a una tantum simulation, the 
perturbative probabilistic distribution built from these input data provides, within an 
approximation with improves with k, the ground state energy E Q as a function of the 
Hamiltonian parameters (coupling and hopping amplitudes). In contrast, in a standard 
Monte Carlo method any different choice of the Hamiltonian parameters requires a 
distinct simulation. Furthermore, the cumulants are easily measured also in the case 
of fermions (bosons in the presence of magnetic fields). The sign (phase) problem that 
occurs in this case remains confined in the expression of the probability distribution and 
can, in principle, be addressed analytically. 

The paper is organized as follows. In section 2 we review the probabilistic approach 
to quantum many-body systems. The equation for the ground-state energy is written 
as an expansion over the cumulants of the potential, hopping and phase multiplicities 
in subsection 2.1 and as an expectation with respect to the probability distribution 
of the same variables in subsection 2.2. In the latter subsection the case of an 
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uncorrelated multinomial distribution is described in detail. In section 3 we introduce 
the probabilistic perturbative scheme bringing together the merits of the multinomial 
probability distribution with the statistical details provided by the cumulant expansion. 
The parameters defining the perturbative probability distribution of the potential, 
hopping and phase multiplicities are explicitly discussed, up to the third order, in 
subsections from 3.1 to 3.5. This is a technical part which could be skipped in a first 
reading. Some considerations on the higher perturbative orders are given in subsection 
3.6. The equation for the evaluation of E resulting from the above perturbative scheme 
is discussed in section 4. Section 5 deals with some numerical results. Concluding 
remarks are drawn in section 6. Two appendices close the paper. In Appendix A 
we review the methods for solving the nonsymmetric algebraic Riccati equation which 
appears in determining the perturbative parameters at the second order. Appendix B 
summarizes the evaluation of the perturbative parameters at the fourth order. 

2. Probabilistic approach to quantum many-body systems 

In this section we review the probabilistic approach to quantum many-body systems 
developed in [9, 10, 11, 12]. Consider a system of particles represented by a Hamiltonian 
operator H acting on a M dimensional space of states labeled by configuration indexes 
n. As an example, think about spinless particles undergoing a simple exclusion dynamics 
in a lattice with configurations given by the site occupations n = (0, 1, 1, 1, 0, 0, 1, . . .). 
In the chosen n-representation, we separate, as usual, H = K + V, where V = diag(i7) 
is called the potential operator and has non vanishing matrix elements 

V n ,n = V n . (1) 

The hopping operator K is defined by the matrix elements 

K n ,n' = —^n,n> Vn,n' , Vn,n> > 0, \^n,n'\ = 0, 1, (2) 

conventionally written in terms of links X n ,n' an d positive strengths r\ n ^ n i. The matrix 
with elements |A„ jn /| forms the so called adjacency matrix and establishes whether two 
configurations n and n' are first neighbours with respect to K or not. The phase change 
(sign in the case of fermions) registered by connecting two first neighbour configurations 
n and n' is registered by Arg A„ jn /. Note that J^ n , \\ n ,n' I is the number of configurations 
which may be connected to the configuration n. 

The evolution of the system from an initial state ipo is earned by solving the 
Schrodinger equation. After a time t the system is found in the state ip{t) which in 
the chosen n-representation has components (we put h — 1) 

(n,V(*)> = ^(n,e^n )(n ,^ ). (3) 

no 

This evolution, as any linear flow, admits an exact probabilistic representation 

[1, 2, 3, 4, 5]. In a one-to-one correspondence with the links, we introduce M 2 

independent Poisson processes {iV£ ,} with rates {p n ,n>}- We recall that for these 
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Figure 1. Random walk with N t jumps in the time interval [0,i): scheme of the 
visited configurations and of the corresponding jump times. 



processes the probability to jump k times in the time interval [t,t + s) is [13] 

Prob(iO - Nl n , = k) = ^^V^, k = 0, 1, 2, . . . . (4) 

We establish that each time a Poisson process N^ n , jumps, the configuration of the 
system changes from n to n' if |A UiT1 '| = 1, otherwise remains n. Arranging the jumps 
according to the times, s\ < s 2 < ■ ■ ■ < s^ t < t < sjv t +i, at which they take place 
in the interval [0,i) we define a random walk in the configuration space of the system 
n — > rix — > n 2 — > ... — > rijv t generated by the above rule from a chosen initial 
configuration no- It is simple to prove that the fundamental matrix elements of the 
evolution operator (n, e~ lHt n ) can be written as an expectation over the above M 2 
independent Poisson processes [14] 

(n,e- i ^n ) = E(5 n , njVt ^^ ), (5) 

Ml = e £".«"W* (^iXkVkPk 1 ^- 1 ^' 3 "- 1 ^ e- iVN ^- SN ^, (6) 

where we put so = and introduced the shorthand 

V k = V n „, k = 0,1,..., N t , (7) 
Afe = A nfc _ 1)nfe , k = l,...,N t , (8) 

Vk=Vn k . 1 ,n k , Pfc=P» k _ 1 ,n fcJ k = 1, . . . , N t . (9) 

The rates {p n ,n'} of the Poisson processes are completely arbitrary, in fact it is easy 
to check that E(dA / l^ /dp riin /) = 0. Here, for simplicity, we take p n ^ n i = p uniform, 
whereas other choices, e.g. p n ,n' — Vn,ri, allow to define optimal Monte Carlo 
numerical algorithms [3, 7]. The above probabilistic representation holds also for non- 
autonomous systems with a time dependent potential V n (t). In this case the time- 
ordered quantum evolution operator Texp(— i J^H(u)du) has matrix elements given 
by (5-6) with Vk — > Jg k+1 Vk{u)du. The representation holds also at imaginary times 
t — > —it with the substitutions A^ — >■ — iA^ and — > — iV^. 

A convenient way to study the properties of the ground state of a particle system 
is to consider its evolution for a long imaginary time. Starting from an arbitrary 
configuration no the system finally relaxes into the ground state, assumed not to be 
orthogonal to no- In this way we can evaluate any ground-state correlation function via 
time asymptotic probabilistic expressions. For instance, the ground state energy E is 
given by 

E = lim -<9 t logV<n,e-^n ) = lim -d t log E(A< ), (10) 
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where -M^ is the imaginary time variant of (6). It is of fundamental importance that 
for the expectation E(A4 t no ) we can find, at large times, an analytical expression. In the 
following we outline how this result is obtained. 

First, we decompose the expectation in a series of canonical expectations over 
random walks with a fixed number of jumps 



E (^o)=E E (- M -o'^ = iV ) 



(11) 



N=0 



and evaluate each term of this series by integrating over all possible jump times. The 
probability to have iV Poisson processes jumping in the interval [0,t) with the k-th 
process jumping in the interval [sk, $k + ds k ) amounts to 

N 

JJ e ~ 'V^-^Vfcdsfc. (12) 

k=l 

With a simple calculation we then obtain 



v 



e (ac N t =N) = Y: w%\t) n A^r, ^ 

ren N k=i 

where Qn = ^v(^-o) is the set of all possible random walks with iV jumps branching 
from n . The contribution of the r-th random walk n — > —>■...—>■ n$ includes 
the factor W$(t) = C'^W^ (z)](t) which is the inverse Laplace transform of 

N 



From equation (13) it is evident that only the random walks with |A 



(r)| 



n k-i 7^ n k'i ioi k — 1, . . . ,N, contribute to the sum. Thus we can rewrite the sum 
over Qn as a probabilistic expectation over these effective random walks. According to 
the choice p n ^ n i = p uniform, the effective random walks correspond to a Markov chain 
with transition matrix P n , n ' = \^n,n'\/Yln" \^n,n"\ [15]. The probability that we must 
associate to the r-th element of the set Qn is, therefore, 



(14) 



1, i.e. 



(r) 

Pn 



A (r) I 



N 



(15) 



Note that p$ = for non effective random walks and J2 r< zQ N P ( N = 1- Multiplying and 



(r) 



dividing by p$, we rewrite the canonical expectations as 

N 



E 

r£Q N 



1 



N 



N 



wri T i r) n A i r) ' ( ig ) 



(r) 



fe=l 



fc=l 



where 



T k = n 



Tifc_i,n|) 



k = l,...,N. 



(17) 



Equation (16) shows that E (Ai* , = ^) is the average of a quantity which does 
not rely on the detailed sequence of the configurations visited during the time t. It 
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depends just on the multiplicities, or numbers of occurrences, of the potential, hopping 
and phase variables, V, T and A, defined by (7), (17) and (8), respectively. For a random 
walk with iV jumps, these multiplicities are explicitly defined as 

N 

N v = Y,$v,v k , VeY, (18) 



fc=0 
N 



N 7 



k=l 



N 



A G J2?, 



(19) 



(20) 



fc=l 



where 2? and J2? are the sets of all possible values assumed by (7), (17) and (8) 
during a random walk with infinitely many jumps. Note that G' J2? as jumps between 
configurations n and n' with \ n , n ' = have zero probability to be realized. Since any 
configuration can be obtained from any other one by a finite number of jumps, i.e. the 
Markov chain we are considering is irreducible, the elements in the sets V, 2? and Jzf 
do not depend on the initial configuration n Q . Let us indicate the multiplicities N v , 
N T and N\ collectively by a vector with as many components as the elements in the set 

/*'' (...Ar-..:....V/ Aa...). (21) 

If we split the set Qn into subsets of random walks with equal values of n, we conclude 
that 



E(M t no ,N t = N) = 



—Nvi 



(t) 11 T Nt ]\ \ N \ 



where 




k-l<' l k 



S.Jr) 



(22) 



(23) 



is the probability to have random walks with multiplicities /j, after jumps from the 
configuration n . Note that Vn{v) = unless satisfies the following three constraints 

J2N V = N+1, ^iV T = iV, J2 N * = N - ( 24 ) 

Vet Te.T \e^f 

As a second step, we observe that when the imaginary time t becomes large the 
full expectation (11) takes exponentially leading contributions from terms with N ~ t. 
Thus, we can replace all previous results with their — > oo asymptotic expressions. Due 
to the ergodicity of the underlying Markov chain, the probability (23) looses memory 
of the initial configuration no- We can evaluate the inverse Laplace transform which 
appears in (22) by a saddle point technique in the complex plane. In the same equation 
we can also substitute the sum over the multiplicities by an integral over fi and avoid 
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to distinguish the normalizations N and N + 1 in (24). As a final result we have the 
following asymptotic logarithm equality 

/ e xot+N(v,u) 
A(Nv) v « (Nv) vm^r (25) 

where we have introduced the frequencies v = fi/N, which have a finite limit for 

iV — > oo, the vectors 

u V (•••• log(x + V) . . . ; . ..log V.. . ; . . . log A. . .), (26) 
•«'' (•••!•'•„ • V) -...:...()...:...()...). (27) 

and the scalar product (a, 6) = J2ae,^' a ab a - The quantity x , which is the real saddle- 
point in the complex contour used to evaluate the Laplace antitransform, is the unique 
solution of the scalar equation 

Z^fv = W x °> - y " < 28 > 

Note that this equation has a regular scaling behavior for t, N — )■ oo with iV ~ t. 

The integral over the frequencies v in equation (25) is easily performed by a saddle- 
point method (actually Laplace's method) if the asymptotic probability density Vn{Nv) 
is known. We cannot hope to evaluate this probability from the microscopic definition 
(23) except for very particular models. For general systems, two different strategies have 
been considered, see [11] and [12]. We can relate Vn{Nv) to proper statistical moments 
of the system under consideration, measure or calculate some of these moments and, 
lastly, obtain partial information about E . Or we can postulate some expressions for 
Vn(Nu) and see which kind of systems are described, exactly or approximately, by this 
guess. 

2.1. Cumulant expansion 

If we rewrite the probability density V^(Nu) in terms of its Fourier transform Vn(q) 
V N {Nu) = (2ir)-^ J dq -p N {q)e-' {q > Nu \ (29) 

we can associate log"Pjv(g) with the cumulants, or connected correlation functions, of 
the multiplicities Nv sampled with respect to the measure Vn(Nv). Indicating with 
{Nv ai . . . Nv ak )ff the component cti . . . of the cumulant of order k, we have the well 
known relation [16] 

logV N (q) = J2y((Nv,iq) k )$ 

k=i 

= Efc! E ••• E (Niy ai ...Nu ak }^q ai ...q ak . (30) 

k=l ' aieJf a k £J$? 

At this point, it is simple to calculate the canonical expectations E (.M* ,N t = N) 
performing the integrals over the 2\J^\ variables v and q by a saddle-point 
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approximation. This approximation is asymptotically exact for N — > oo. The last 
step is to evaluate E by resumming the series (11). By virtue of the limit t — > oo to 
be taken at the end, we still have an exact result if we replace the series over N by an 
integral and estimate this integral at its saddle-point. There is one important point to 
observe. Independently of their order k, the cumulants (Nu ai . . . Nv ak )^ diverge as N 
for N — > oo. Thus we introduce the asymptotic rescaled cumulants 

^l.,a k -^{N ai ...N ak )^, (31) 

£( fc ) in a compact notation. Existence and finiteness of these limits are ensured by the 
finite correlation length N c which characterizes the correlations functions, connected or 
not, of the the multiplicities [11]. We conclude that E is the unique solution of the 
scalar equation 

E jfci £ • • • E E £L« h ««i(^) • • • u ak (E ) = 0, E < V min , (32) 

k=l ai£,34? OLfcdffl 

where 

u T (E ) = (...- log(-£ + V) . . . ; . . .logT. . . ; . . . log A. . .). (33) 

Equation (32) is exact and the uniqueness of its solution is ensured by the constraint 
E < V ni i n which stems from the Laplace transform causality condition (28). 

To find the exact ground-state energy E from equation (32) we have to know the 
cumulants T,^ at any order k. Of course, for a general system this is not conceivable. 
However, up to some small order k and even for systems of relatively large size, the 
cumulants can be measured by reliable statistical simulations [11]. With these input 
data, we can truncate the series in (32) at some order and obtain an approximation 
to Eq. Independently of the truncation order, there is an important feature to note. 
Suppose that we change the Hamiltonian H leaving unaltered the adjacency matrix 
| A™,™' | and the number of elements in the sets Y , 2? and ££ . The asymptotic rescaled 
cumulants are unaffected by this change and the only modifications to equation (32) 
are encoded analytically by u(E ). Therefore, the same input data S^, k — 1, . . . , /c max , 
can be used to find the ground-state energy of parametric Hamiltonians as a function 
of their parameters. The simplest example is to consider H(^y) = K + and evaluate 
the function Eq(^). 

At the lowest truncation order A; max = 1, equation (32) reads 

J2 4 1} M-£ + v) = J2 s ? io g T + J2 log A > 

vet Tasr \eJ? 

E < V min . (34) 

In general this equation must be solved numerically. Conversely, for V = we have 
Y = {0} and Sy^ = 1 which allow us to find the analytical solution 

4 0) = -[ !! ! ( 1 1 \ ]. (35) 
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This represents the lowest order approximation to the ground-state energy of the 
hopping operator K. 

The cumulant expansion described above has been applied to study Hubbard models 
in a two dimensional lattice [11]. By including cumulants up to order 4, the results 
compare rather well with the exact ground-state energy (determined numerically in 
other ways) at least in the case of hard core bosons, i.e. ££ = {1}, and for interaction 
energies not too large with respect to the hopping term. The results, however, are 
disappointing for large interaction energies, in fact in this limit E diverges. Moreover, 
in the case of fermions when J*f = { — 1, 1} pointless complex solutions can be found 
for E as shown, for instance, for V = and at the lowest order by equation (35). 
These problems stem form the fact that truncating at some order k max equation (32) is 
equivalent to consider an artificial probability density Vn(Nv) with zero cumulants at 
any order k > fc max . 



2.2. Multinomial probability density 

In view of the definition (31) of the asymptotic rescaled cumulants, equation (32) can 
be compactly rewritten as 

oo 

JToo N £ M ((A ^' U W)*>" = ' E ° ^ (36) 

00 k=l 

The series sums up to logT^— iu(E )) hence we can state that E is the unique solution 
of the exact equation 



lim log / d(Nu) V N (Nv) e^-^)) = 0, E < V min , 



(37) 



with u(E ) given by (33). Equation (37) makes crystal clear that the knowledge of E 
stems from that of Vn(Nv). 

The random walks in the configuration space definitely induce correlations among 
the multiplicities /j, due to the dependence of the potential, hopping and phase variables 
on the visited configurations. Could we neglect these correlations, for iV sufficiently 
large the multiplicities of each set V, 2? and ££ would be equivalent to multinomial 
trials processes with success probabilities 

1 N 

Pv = £ m oo iUl £ ^ = V e r, (38) 



N^oo N + _ 

1 N 



^ = &,^E^=4 1) , Te^r, (39) 



k=i 

N 



PA= J fc,^E^ = 4 1) . ^y, (40) 
k=i 

Note that we made use of the ergodic properties of the underlying Markov chain. In 



A perturbative probabilistic approach to quantum many-body systems 



11 



that case, the probability Vn (/•*) would be a product of multinomial distributions 

Ny N T N x 

where the multiplicities N V) N T and N\ are integers which satisfy the constraints (24). 

Equation (41) greatly simplifies for N large. By using Stirling's approximation 
for the factorials and explicitly taking into account the constraints (24), we have the 
following asymptotic equality for the associated probability density. 



V N {Nu) ~ exp{Ncu(is)] 5 V Nv v ~ N 



x 5 



| Nu T - N J 5 ( E ^a - N J , (42) 
\Te5* / \Ae^ / 



where 



(43) 

and p T = (. . . py px ...;... Pa • • •) is a vector collecting the success probabilities 

(38), (39) and (40). To find E it remains to calculate the integral of equation (37) over 
the frequencies v. The integration can be performed by steepest descent as detailed in 
the following section. The result is that E is the solution of 

^ -E + V (Y,Te? PT T ) [ExeSf PaA) 

For V — 0, i.e. ^ = {0} and pv=o = 1, equation (44) can be solved analytically and we 
obtain the value 

4 0) = -(5>t) (j2p x x) (45) 
\re^ / \Ae^ / 

for the ground-state energy of the hopping operator K. For V ^ equation (44) is 
straightforwardly solved numerically by bisection method. 

The uncorrelated multinomial probability density considered here has the great 
advantage that cumulants of any order are included in the determination of E . This 
eliminates the artifacts obtained by truncating the cumulant expansion (32), namely 
the wrong behavior of E at large interaction energies and its complex value in the case 
of fermions. In fact, for large interaction strength equation (44) admits the asymptotic 
finite solution 

E = V min + pv^E™ . (46) 

Moreover, for intermediate interaction strengths the solution of equation (44) varies 
monotonously between the two limits (45) and (46), as expected. In the case of fermions, 
the V = solution (45) is real and negative, as must be. In fact, among the links 
Aji, n ' 7^ those with a positive sign are the majority so that 

5>A = E£i " > 0. (47) 

Ae.Sf 
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The main drawback of the present approach has already been mentioned. The 
correlations among the multiplicities p, are neglected so that the comparison of E , 
solution of equation (44), with the ground-state energy of a system of particles can 
be quantitatively poor. Whenever these correlations are absent, as in the case of the 
uniformly fully connected models and of the random potential systems considered in 
[12], the approach provides exact results together with the possibility to study the 
appearance of quantum phase transitions in the thermodynamic limit. 



3. Multinomial perturbative scheme 

Here we propose a probabilistic perturbative scheme with the aim of merging the merits 
of the multinomial probability density of section 2.2 with the statistical details provided 
by the cumulant expansion of section 2.1. Basically the idea is as follows. We consider 
an asymptotic probability density which has the same structure of a multinomial density 
but with parameters p a , a & JF, that are functions of the frequencies v>, 



co (u) = lim 1 \o g V N (Nu) = V v a log 



(48) 



and write these functions as power series of the form 

oo 1 

*■ <"> = E h E ■ ■ • E & (>* - p { £) • • • (•*. - pg.') • (*» 

n=0 ' /3i6JT /3 n eJ? 

The scalars p^l. ,a k , with Oil E for i — 1, . . . , k, in a compact notation p( k \ are the 
perturbative parameters at order k. Note that we identify the perturbative order k with 

the rank of the tensor p^ not the index n of the series (49), i.e. k = n + 1 = 1,2, 

At the lowest perturbative order k = 1 we have p a (u) = pa \ a G Jtf, constant as in 
the strict multinomial case. 

We determine the perturbative parameters as follows. First, we assume that 
for k > 2 they are symmetric under the exchange of any two of their indices eti, . . . , 

P%... aj ..= P%...a,..- (50) 

Second, in order to maintain as much as possible the structure of a multinomial, we ask 
that, for any value of v, the functions p a (u) are normalized in each set Y, i.e. 

= 1, ^ = r,^,j2?. (51) 

Third, we require that the asymptotic rescaled cumulants of order k evaluated from 
(48-49), hereafter indicated by ({u ai . . . v a/ J), coincide with those effectively owned by 
the system, 

K 1 ...i/« fc >>(p (1 \p (2 \...,P ( * ) ) = E£U- (52) 
In this expression we have anticipated that ((z/ ai . . . v a ^t) depends only on the parameters 
with j < k. This property implies that we can first find p^ by solving the system 
of \Jf\ equations 

K.Mp^^eW, ai eJt?, (53) 
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next find p^ by solving the system of \Jj? | 2 equations 

iv ai v a2 W l \p {2) ) = SgU. (54) 
next find p*- 3 ^ by solving the system of |^| 3 equations 

((iy ai ^a 3 ))(p {1 \p {2 \p {3) ) = SgUa' «1,«2,« 3 e , (55) 

and so on up to a chosen maximum order /c max which corresponds to truncate the series 
(49) at the term n = k max — 1 included. 

To be precise, the probability density (48) is well defined only if each p a (u), a G Jff, 
is a non negative function of the vector v with components v a > varying in the unit 
simplex ^ agrf = 1 for srf = Y, ST . Actually, we will content ourselves that the 
positivity condition for p{y) is satisfied in a small region around the frequency v which 
gives an exponentially leading contribution to the integral (37). A priori we don't know 
whether this condition can be met by introducing the perturbative parameters p^ as 
described above. We will then proceed heuristically. Whenever an effective solution of 
equation (37) can be found by using the probability density (48-49), that will be the 
signal that the proposed perturbative scheme is meaningful. 

In the following subsections, we will first evaluate, up to the third order, the 
asymptotic rescaled cumulants associated to the probability density (48-49), then we 
will discuss the solution of the systems of equations (53), (54) and (55). Finally, we will 
make some comments on the higher perturbative orders. 



3.1. Evaluation of the asymptotic rescaled cumulants 

The definition of the asymptotic rescaled cumulants of order k — 1, 2, . . . is 

({u ai ...u ak )) = Jim (Nu ai . . . Nis ak ) { $ , a u . . . , a k e Jf, (56) 

where (Nu ai . . .Ni/ Qk )ff are the cumulants, or connected correlation functions, of the 
multiplicities N a = Nu a sampled with respect to the iV-jumps probability density 
Vn(Nv). In turn, as standard, the cumulants are obtained from the generating function 
associated to Vn{Nu) 



{Nv ai ...Nv ak y N 



(c) _ d k log Z N (J) 



dJ ai . . . dJ ak 



(57) 



j=o 



Z N {J) = J d(Nv)V N (Nv)e( J > N »l (58) 

Assuming that for — > oo the logarithm of Vn(Nu) is given by (48-49), for N large 
we have, up to an inessential constant, 



Z N (J) = I d(Nv)e N "^ + ( J ' N ^ J] 5 ( Nu a - N J , 



(59) 



where we have explicitly taken into account the fact that the multiplicities of each set 
y, 3?, =Sf must sum to N. Using the Fourier integral representation of the Dirac 8, we 
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rewrite the generating function, up to a constant, as 

Zn{J)= [ J] dk^ H du a e N ^ J \ (60) 

where 

<P(u, k, J) — "* flog ( — ) + J « 

+ Yl (61) 

In the above expressions we used the compact notation k = {ky , k & , k%) and J = 
(... J v J T J x . . .) with V e f , T G J and A 6 ^. Evaluating the 

integrals in equation (60) by Laplace's method, we get the iV-large asymptotic logarithm 
equality 

Z N (J) ~ e N ^ J \ sp ( J) = (z/ sp (J), fc sp (J), J) , (62) 

where (z^ 15 , fc sp ) is the saddle point (actually, the maximum point) of <fi solution of the 
following system of equations 

0, aeif, (63) 



Jj^ = 0, */ = r,&,J?. (64) 
By differentiating equation (61), the saddle-point equations can be written as 

(,.sp \ D ' (l/ sp ) 

) = j tt + y ^ a 7 ; ^ p +i^-i, (65) 

E^ P = 1 ' ^ = r,3T,^, (66) 
where we have defined 



l>S+El>S,K- *?>)+" (67) 



and 





a G f, 






-2\ 


a e JSf. 



(68) 

Note that, due to the properties assumed for the perturbative parameters, p' aj3 (u) is 
symmetric under the exchange of the indices a, (3. From equation (65) we find 

/ / sp\ 

= Pa („<*) e Ja+ ^~ P ,(- sp ) +lk ^ \ aG/, (69) 
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which inserted into equation (66) provides 

iA- sp -l 1 

e lfi w 1 — 



15 



p' «(- sp ) 



(70) 



By using equation (70) we conclude that the saddle-point frequencies (69) are the 
solution of the system of nonlinear equations 



a 



Pa \y° v ) C 



sp\ ^ J "^2^/3eM- P p(v s P) v p 



(71) 



Note that both i/ sp and fc sp are functions of the source J. It follows that the function 
4> (v, k, J) evaluated at the saddle point (u sp ( J), k sp (J)) is 



K (J) = E 



aeJf 



log 



Pa [y 

sp 
Va 



+ Ja 



= E iog(E^^ sp ) eJa+E ' 



^(- SP ) ..sp N 



We are now ready to evaluate the asymptotic rescaled cumulants associated to the 
probability density (48-49) by using the formula 

<9 fe 0s P (J) 



{(u ai ...V, 



0J ai . . . 0J ak 



(73) 



J=0 



3.2. Derivatives of <p sp (J) 

In this subsection we supply the derivatives of <p sp (J) with respect to the source J up 
to the third order. Note that sp depends on J explicitly and through the saddle-point 
frequencies. 

The first order derivative of 4> SV (J) with respect to the source component J ai , 
ai E Jt?, is 



<90 sp (J) 



dJ. 



ai 



E E"? 



1 dp a (u sp ) 



p a {v sp ) dJ t 



+ 5, 



ai 



-EE 



dv S J> P' a0 (^ SP ) 



sp 



dJ ai Pn (v sp ) 



Since 



dp a {u sp ) 



dJ, 



E ^ (^ sp ) 



dJ, 



and p' a8 is symmetric, the previous expression reduces to 



d(p sp (J) 



v. 



sp 

Ql ' 



(74) 



(75) 
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The second order derivative of 4> sp (J) with respect to the source components J ai 
and J a2 , ai,a 2 G Jf, is 

d 2 sp (J) _ dv^ 
8 7 87 8 7 

= U SP S -Z/ SP Z/ P Y 

1 8 Pa (^) ^ 8 (p' aP {v sv 



x 



Pa(" sp ) dJ a2 + J^ f dJ a2 {pp(v^ l/ 



where 



= n otherwise . ( 76 ) 



Recalling the definition (67) and using equation (74), we find 
d 2< kp( J ) = x _ ^sp s P 

(7 J Ql OJ a2 

<9z/ p 

+ EE - « P w) (^ sp ) (77) 

where we have introduced the symmetric matrix 

and defined 

8 2 p a {iy) 



= pS, + E !>$>.(•*- + - (79) 

Note that, due to the properties assumed for the perturbative parameters, p" a ^(y) is 
symmetric under the exchange of any two of its indices a, (3, 7. 

The third order derivative of </> Bp (J) with respect to the source components J ai , J a2 
and J as , ai,Oi2,Oi3 G Jff, is 

(97 (97 (97 (97 87 

8v sp fdv sp 8u sp \ 

= ai 5 - { ai U SP + Z/ SP " 2 I Y 

(m ( sp ) + 9M ^^ SP ) 9 ^ P> \ ( 80 ) 
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with 

dM aP (z/ sp ) 



dJ, 



a 3 



E 



p'Lm (^ sp ) , p'Lm (^ sp ) + p'Lm (^ sp ) 



Pa (^ SP ) Pfi (^ sp ) p 7 (^ sp ) 

P'afS K 7 (^ SP ) K/3 (^ SP ) P'M (^ SP ) K 7 (^ SP ) P' lP (" S " 



Pa (V SP Y 



pv iy sp Y 



p~< (" sp Y 



PaW^) P'^S^P's,^) 

Ps (^ sp ) p 5 iy^f 

PaS, (^ SP ) P'sp (^ SP ) P'SM (^ SP ) PL (^ SP ) 



ps {v s *y 



Ps (^ sp )' 



+ 2 

where we have defined 



P'aS (^ SP ) ^ (^ SP ) 



PS (^ sp ) 



sp 



d]y sp 



dJ, 



(81) 



Q3 



= P ( Xs + J] _ P^) + • • • 



(82) 



3.3. Equations for the perturbative parameters: first order 



The perturbative parameters p^ are determined by the system of equations (53) which, 
according to (73) and (75), reads 

(83) 



v. 



sp 

ai \J= 



=0 = s £> «i e ^- 

Using equation (71), we explicitly have 

pL i/3 (s (1) > 



p ai (S«)e 



E 



(i) 



^•^ P(3 (s(i)) /3 



,(=(!)). 



El 1 ) 



ai ) 



(84) 



A solution of this system of nonlinear equations is 

P# = Eg, «i e (85) 

To prove it, we start to observe that if p^ = E^ then for any a, (3 E Jt? 

Pa(W)= P £\ (86) 

and 

P'a,^) = P { al (87) 

The position p^ = E 1 ^ has another consequence which is pivotal also in determining 
the perturbative parameteres of higher order. In fact, due to the constraints (51), the 
analogous normalization conditions valid for E^ and the sysmmetry properties (50), 
we have 

k>2, 1 < i < k, £/ = y,&,Sf. (88) 



V v (k) = o 
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In the present case, the above sum rules allow to write pes? Pap = 0> wmcn > together 
with (86) and (87), permits to reduce the l.h.s. of equation (84) to . 

Equation (85) is not the unique solution of the system (84). However, besides 
being the natural solution for which the perturbed probability density (48-49) reduces, 
for k max = 1, to the uncorrelated multinomial (43), it allows us to establish the sum 
rules (88). We shall show that these, in turn, entail the asymptotic rescaled cumulants 
((u ai . . . v ak )) of arbitrary order k to depend only by the parameters p^ with j < k. In 
this way, the parameters p^\p^ 2 \p^\ ... we find do not depend on the order k max at 
which we decide to truncate the series (49). 

3.4- Equations for the perturbative parameters: second order 

The perturbative parameters p^ are determined by the system of equations (54) which, 
according to (73) and (77), reads 



(2) (2) / (3) (2) (2)\ 

P»P + PaP_ + / PaP 1 _ PaiP 1 p \ ^(1) 



,(1) ' y(l) ' \ y(l) (1)2 J 

J " ^P T 6JT V^T ^7 / 



S (2) = V(2) 
Pa.2 ct\ct2' 



E (2,0) + S (2,0) 

In writing this expression we have used the results (83), (85), (86) and (87) of the 
previous section as well as the fact that, since p^ = S^ 1 ', for any a, (3, 7 G Jt? it is 

P'Lp,^ {l) )-pfpr (89) 
We have also defined 

SS 0) = ^5 a p " ^a^Xap, arfeJT, (90) 

E^ 2,0 ) in a compact notation, which is the second cumulant of an uncorrelated 
multinomial probability density with parameters E^. According to equation 
we have X^e^r Pap-y = so that the above system of equations becomes 



,(2,0) 
J aia 



(2) (2) (2) (2) 

Pap PaP _ ST^ P a ~fPyP 

+ ~ 2^ y(l) 

Z-'fj ieJ #> ^7 



y(2) _ y(2) _ y(2,0) / q-i \ 



For 0:1,0:2 G Jf, this is a nonlinear system of m 2 , m = \Jf?\, equations. On the other 
hand, p( 2 \ recall that it satisfies the sum rules (88), has m 2 , m = — 1) + (|^| — 
1) + (|J2?| — 1), independent components. Therefore, the system of equations (91) is 
overdetermined and we have to lower its dimension to find p( 2 \ Note that also the 
m x m matrices E^ 2 - 1 and E^ 2 ' -* are singular and their rank is m 2 . 

Let us introduce the m-dimensional index space J4? = J$? \ {V*,T*, A*}, where K, 
T* and A* are three arbitrarily chosen elements of the sets Y, and J2? , respectively. 
Let be the vector with components p^ 1 = Pa\ ot G and p^ the matrix 
with components p^ = p^ 2 j, a, (3 G J&. Now consider the rh 2 equations (91) with 

01,02 G ffl ■ Observing that E Qg ^ = ^ tt e,f e{v*,T*,A*} an d usm g the sum rules 
(88), we rewrite the first term in the l.h.s. of equation (91) as 

J 2 ) / V (2,0) y (2,0) \ 

v (1) Zj /3«2 ~ 2^ I v (l) " v (l) I P aP [^fa ^P*a 2 ) ' 

ae^pGJT aejrpej? \ ^ a / 
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where a*, f3* are the components eliminated from the sets &>? a , &/p, i.e. 





[ K, 






a* = | 




a G ^, 


(92) 




[ A*, 


a G J2f. 





Similarly, the second term in the l.h.s. of (91) becomes 

(2) /y( 2 ) y( 2 ) 

y(2,0) P<*P y(2) /v,(2,0) _ W2,0)w(2) f _ f^ag 

Z^ Z^ «i« v (l) /?"2 _ Z^ Z^ I "l« ^a ia J Paf3 I (i) (i) 

whereas the third term gives 

(2) (2) 

E>yp v (2,0) gg2^. y( 2 ) 
Z^ Z^ ^"1" v (l) 

ae,je pejtr jejtr ^7 

= £EIE M? - f & + ») #S (4 2 i - 41) ■ 

aej? pejr -yejr sej? V 7 7 * / 

We conclude that p( 2 ) is the solution of the nonlinear matrix equation 

S (2 ' 0) p(2)E( 2 ) + S( 2 'V 2 )S (2) - Y.^p^Tp^Y.^ = S^ 2 ) - S( 2 '°). (93) 

where £( 2 ) and £( 2 '°) are the reduced versions of the matrices £( 2 ) and £( 2 '°) and we 
have introduced the matrices Y,( 2 \ £ , £( 2, °), £ ' and T with components a, (3 E 
given by 



^(2) v (2) v (2) 


(94) 


v (2) v (2) 
VT(2) ^a*0 

^ y(l) V (l) ' 


(95) 


^(2,0) _ v (2,0) y (2,0) 


(96) 


v (2,0) v (2,0) 
vr(2,0) ^a/3„ 

a ^ ~ V (l) V (l) ' 
^P ^P* 


(97) 


r <W Xa/3 


(98) 



Since 2V 2 ) and £( 2 >°) are nonsingular, equation (93) can be rewritten as 
S(2-o)- 1 s (2 'V 2 ) +p(2)s (2) vM-'-p^TpW 

= S( 2 <°) ~* (t^ - S( 2 '°)) S( 2 ) _1 . (99) 

This is a nonsymmetric algebraic Riccati equation (NARE), which can be solved 
numerically by matrix factorization (Schiir method) [17] or iterative methods [18]. For 
details we refer to Appendix A. Once p^ has been found, the complete set of second 
order perturbative parameters p^ is obtained using the sum rules (88) for k = 2. 
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3.5. Equations for the perturbative parameters: third order 

The perturbative parameters are determined by the system of equations (55) which, 
according to (73) and (80), reads 

/_(2) (2) / (3) „(2)_(2)\ \ 

V(3,0) , y(3,0) f Pap_ P*P_ / PaPy _ P^hl \ yW I v( 2 ) 

«i«2«3 ^ Z^ Z^ «i««3 I v (l) ^ v (l) ^ I v (l) (1)2 I ^7 I ^/3« 2 

aeJf/SeJT V^" ^fi 7 6JT V ^7 L 7 / / 

/_(2) (2) / (3) „(2) (2)\ \ 

^(2,0) [ Pa0_ Paj3_ Sp I Pap-y _ PcnP^p \ yj(l) j yj(3) 



+ Z^ 1 v (l) + v (l) t Z^ I v (l) ~ (1)2 I ^7 I ~/Ja 2 a 3 

aeJT/3e^ V^" ^ 7€JT \ Zj 7 ^7 / / 



+ E E s 2 E 



(3) (3) (3) (2) (2) (2) (2) (2) (2) 

Pa^-y Pa/3-y Pa/3y _ Pp a P a l _ PafiPfiy _ P a lP~f/3 

WWW "v^" y^~ 

^7 2_, a 2^ 2^ 7 



+ E 



(4) (3) (2) (3) (2) (3) (2) 

PgfiyS _ PafiSPs-y _ PaSyPsp _ Psp-yPSa 

w w ~W 2 w~ 



S&JT \ ^5 ^ s ^ s ^ s 

(2) (2) (2) 
} Pa8P(5SPy8 \ v (l) 



2 J "° J pu : < u n 

4 



ZflMs = S £Ua- ( 10 °) 



In writing this expression we have used the results (83), (85), (86), (87) and (89) of the 
previous sections as well as the fact that, since p^ = E^, for any a, 7, S G Jt? it is 

P'^s(^ 1) )=P%- (101) 
We have also defined 

E$? = ^ - + Xa p, (102) 

E^ 3,0 ^ in a compact notation, which is, formally, the third cumulant of an uncorrected 
multinomial probability density with the first two cumulants equal to E^ and E^ 2 \ 
respectively. According to equation (88), in the first two lines of (100) we have 
^2-yejf Pa/3-y = an d i n the fourth line ^2seJt? Pap-yS = 0> so ^hat equation (100) becomes 

V V V fA( 2 <°)£ (2) E( 2 ) + £( 2 '°)A (2) S( 2 ) + £( 2 >°)£ (2) A( 2 ) W 3) 

a£,je peje -yeje 

where we have introduced the matrices A^ 2 ^ and A^ 2,0 ) with components ct,/3 e 

' ' ' ' 1 ' (104) 



= f ^ ~ E « ] ^T, (105) 




and the tensor A with components cci, «2, «3 G 
A = S (3) - S (3 - 0) 

t - i fi"2a3 ^010203 ^010203 
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- V V s (3 ' 0) 



(2) (2) (2) (2)\ 

Pa/3 Paf3 _ ST^ P^P-yP \ 

+ yW 2^ 

^P 



y(l) 



yW 

■yeJT ^7 



) 



y(2) 



V(2,0) 



Pal Pal 
yW + yW 



(2) (2)\ 

E P^PfP \ y(3) 

v (l) J ^P^aa 
7 6JT ^7 



+ EEE S S 



(2) (2) 
PWpPo.1 



yV 



(2) (2) 
PapPpJ 

yW 2 
^P 



+ 



(2) (2) 
PajPjP 



- 2 E 



(2) (2) (2)' 
PaSPpsPyS 



S (2) V(2) 



(106) 



Note that, once p^ and p^ have been determined, this tensor can be considered known. 

Equation (103) is, for ai, 012,013 G Jif, an overdetermined linear system of m 3 
equations in the unknown p^ which, according to the sum rules (88), has m 3 
independent components. To find p^ we proceed as in the previous section. Let p^ be 
the tensor with components p a p y = P a p-yi a iftil G J^, and consider the m 3 equations 

(103) with a 1 ,a 2 ,a 3 G Observing that J2 a eJf = Saejr + Z)ae{v„T.,A,} and usin § 
the sum rules (88), we rewrite equation (103) as 



V V V U(2,0)y(2) y(2) ,y(2,0)7(2) ^(2) , f (2,0)y(2) T(2) \ -(J 



(3) 

Pi 



aejt? peJf -yeje 

^aia 2 «3; (^07) 

where we have introduced the matrices A^ 2 ) and A^ 2 ' ) with components a, f3 G 

AS = AS-AS„ (108) 
A^ 0) 



A^-A^, 



(109) 



whereas the matrices £( 2 ) and E^ 2 ' ) are defined by (94) and (96). 

The system of equations (107) with oti, a 2 , a 3 G is a linear tensorial equation in 
the unknown p( 3 \ It can be solved by vectorization, i.e. by introducing a bijective map 
between the set ,3& 3 and the integers {1, 2, ... , m 3 }. Let n(a) : J$? h- >■ {1,2,..., m} be 
some ordering of the elements a G Jt? and : {1,2, .. . , m) h- >■ its inverse. We 
define the integer map i(a, /3, 7) : ^ 3 >->■ {1,2,..., m 3 } by 

i(a, /3, 7) = (n(a) — l)m 2 + (n((3) — l)m + 71(7), 

so that the inverse triplet (a(i), f3(i), is given by 

i-0+ 



(110) 



a(i) 

7(0 



n 



n 



n 



+ 1 



i - (n(a(i)) - l)m 2 - 0+ 



m 



+ 1 



1 (i - (n(a(i)) - l)m 2 - (n(/3(i)) - l)m) , 



(111) 

(112) 
(113) 
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where + is a positive infinitesimal. Instead of equation (107) with a±,a2,as G Jtf, we 
thus obtain the equivalent system 

rh 3 

^Qi, pf } = A>, t = l,2,...,m 3 , (114) 

3=1 

where we have defined 

n - a( 2 '°) V( 2 ) V( 2 ) i v(2.°) A ( 2 ) V(2) , y(2,0) Pi(2) T(2) n . _v 

Vij - 'V0«(j') WMO 7CihM + ^a(i)a(j) PU)P(i) lUhtt) + "(iM^WW 7007(0' ^ ii0 ^ 

as well as p) = Pa(j-),80'b0') and A; = A QW/ 3 W7(i) . Equation (114) is a linear matrix 
equation which can be solved by standard methods, e.g. LU-factorization [19]. 

Once p^ has been found, the complete set of third order perturbative parameters 
p^ is obtained using the sum rules (88) for k — 3. 

3.6. Some considerations on higher orders 

In the previous subsections we have evaluated the perturbative parameters p^ for the 
first three perturbative orders k = 1,2,3. In this section we will show that for any 
k > 3 the equations which determine p^ are i) linear tensorial equations ii) depending 
only on the parameters p^ with j < k. This ensures that the perturbative parameters 
are, in principle, solvable at all orders and that their value is independent of the integer 
kma, x — 1 at which we decide to truncate the series (49). The evaluation of p^ is outlined 
in Appendix B. 

To demonstrate the properties i) and ii), first of all let us point up why they 
hold in the analyzed case k — 3. The tensor pa}a 2 a z is determined by the system of 
equations (55) the structure of which is established, see equation (73), by the third 
order derivatives d 3 4> sp (J)/dJ ai dJ a2 dJ a3 evaluated at J = 0. According to equation 
(80), these derivatives contain rational combinations of the functions p(u) and their 
derivatives p'(v), p"{y) and p"'(y\ a compact notation to indicate respectively the 
components (49), (67), (79) and (82), evaluated at v — u sp . When we set J = 0, since 
v sp \ J=0 = and p {l) = £ (1) , we have 

P(» sp )\j =0 =P W , (H6) 
P(» sp )\j = o=P (2 \ (H7) 
P"(» SP )\ J=0 =P (3 \ (H8) 
p"V P )| J=0 =P {4) , (H9) 

and so on. Notice that in the third order derivatives (80) there is only one term which 

contains p'", namely 

When this term is evaluated at J = 0, the factors vf\ J=0 and p s (t /Sp )| J=0 cancels each 
other out so that, according to the sum rules (88), we have 
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As a consequence, the system of equations which determines p^ does not contain 
p( 4 ). The property is immediately extended to higher orders. The system 
of equations which determines paia 2 a 3 a 4 depends on the fourth order derivatives 
d 4 (j) S p(J)/dJ ai dJ a2 dJ a3 dJ a4 and these contain p, p', p" ', p'" and p"". The fourth order 
derivative p"" may appear only in the term 

obtained differentiating the factor p'"^ s (u sp ) of (120) with respect to J a4 . When 
evaluated at J = 0, equation (122) vanishes and we find that no p^ terms are contained 
in the system of equations for p^. Iterating, we conclude that the property ii) is valid 
at any higher order. 

Now we focus on the property i). The system of equations which determines 
p( 3 ) is linear in the unknown tensor simply because the third order derivatives 
d 3 4> sp (J)/dJ ai dJ a2 dJ a3 are linear in p" ', see equation (80). By further differentiating 
(80) with respect to J Q4 there is no possibility to generate a term nonlinear in p'", 
e.g. quadratic. In fact, this would amount to have in d 3 (f) sp (J) /dJ ai dJ a2 dJ a3 a term 
containing the product or the ratio of the components of p'" and p" . However, we have 
seen that the only term of (80) containing p'" is given by (120). We conclude that the 
system of equations which determines p^ is linear and, iterating, the same holds at any 
higher order. 



4. Evaluation of the ground-state energy 

We have already seen that the ground state energy E is the unique solution of the 
scalar equation 

lim -J- log In(E ) = 0, E < V miQ , (123) 



where 



and 



I N (E ) = J d(Nu)V N (Niy)e^ u ^\ (124) 



u T (E ) = (...- log(-£ + V) . . . ; . . .logT. . . ; . . . log A . . .). (125) 

In this section we want to determine Eq when the probability density Vn(Nv) is given 
by the multinomial perturbative scheme (48-49). 

For iV large we have already calculated the integral (124). In fact, I N (E ) coincides 
with the generating function Zn(J) studied in section 3.1 provided we choose the source 
J = u(E ). Therefore, equation (123) is equivalent to 

<P sp (u(E )) = 0, E <V min , (126) 

where (f) sp (u(E )) is given by (72) with J = u(E ). Of course, also the saddle-point 
frequencies (71) which appear in (72) must be evaluated with the same choice of J. We 
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conclude that the ground-state energy E Q is obtained, together with the saddle-point 
frequencies v sp , as the solution of the following system of nonlinear coupled equations 

v v *V" sp > L .sp r sp 

h ZI ^+ V= (E t ^pA^)t) (E«paWA)' Eo - Fmin ' (127) 

z/p = Pa ^ ae/, (128) 

where we have defined 

p Q H=p Q (^)e E ^W^, (129) 
and we recall that 

* <"> = E ^ E ■ • • E At'k - 4?) ■ ■ ■ K - =S) . 

n=0 /3i6jT ,3 n €Jr 
Pa/3 ( |/ ) = X] ( n _ 1)! • • • Pa'Si-Tn-i 

xK-EW)...K.,-EW l ). 

Note that the above system of equations is valid at any perturbative order. More 
precisely, for different choices of k mSLX only the functions p a (u) and p' aj3 (u) are to be 
modified whereas the structure of the equations (127-128) remain unchanged. 

At the lowest perturbative order k max = 1, we have p(u) = E« p'(v) = 
and p{y) = EM. It follows that equations (127) and (128) can be solved separately. 
Equation (127) reduces to (44), i.e. E is the solution of 

S M i 

™ - Eo + v (£t 6 * 4^) (Ea G ^ 

This equation can be straightforwardly solved by bisection method. Once E is found, 
the saddle-point frequencies are given by 

^(-Ep + V)- 1 



4 P = - v ^° l v g r, (i3i) 



y(l) T 

4 P = T m , T G (132) 

4 = m , A G (133) 

Z^A'eif A' A 

At higher perturbative orders k max > 1, once the perturbative parameters p( fc ) 
for k = 1, . . . , k max have been determined from the corresponding exact cumulants 
EM, . . . ; $]( fcmax ), the system of equations (127) and (128) can be solved numerically by 
a globally convergent multidimensional Newton-Raphson method [19]. Of course, we 
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have to conjecture some initial value of (E Q , v sp ) which is not too far from the solution. 
Actually, this may not represent a problem in the spirit of our perturbative approach 
according to which the uncorrelated multinomial probability density can roughly capture 
the features of the system considered. Therefore, we propose to use the solution of (130) 
and the values (131-133) as an initial guess for (E ,v sp ). 

5. Numerical results for the Hubbard model 

In this section we apply the multinomial perturbative scheme to study some example 
systems. We will focus on the Hubbard model with spin Vi hard-core bosons. The case 
of fermions will be considered in a paper dedicated to the sign problem. 

The Hubbard model is one of the simplest models displaying the real word features 
of a many-body system. It plays essentially the same role in the problem of the electron 
correlations as the Ising model in the problem of spin correlations. The model describes 
interacting particles, bosons or fermions, in a lattice, and the corresponding Hamiltonian 
operator consists of two terms: a kinetic term allowing for hopping of particles among the 
sites of the lattice and a potential term consisting of an on-site interaction. The model 
was originally proposed by Hubbard [20] to describe electrons in solids and has since 
been the focus of particular interest as a model for high-temperature superconductivity. 
Recently a large interest has been devoted also to the properties of the 2D Hubbard 
model on the honeycomb lattice, as a basic model for the description of graphene [21]. 

Let us consider the first-neighbour uniform (FNU) Hubbard model defined by the 
Hamiltonian 

k = E E ( £ ts> + sW) + t E a It a *t a lA' ( 134 ) 

(i,j)er<r=n *eA 

where A C Z rf is a finite rf-dimensional lattice with |A| ordered sites, T = i < 

j G A and i,j are first neighbour}, and c ia a commuting destruction operator at site 
i G A and spin index a =|! with the property cf a = (hard-core boson destruction 
operator). Hopping strengths are uniform through the lattice and described by the 
parameter r\ > 0. Also the on-site interactions are independent of the site index and 
their value is fixed by the parameter 7 > 0. The system is described in terms of Fock 
states labelled by the configurations n = (n^n^, . . . , 7iiA|t> ^IAll) with n ia = 0, 1. In 
this n-representation, the on-site Hubbard interaction is the diagonal potential operator 
V with matrix elements 

V n , n = V n = 7 ^ n^n ih (135) 

ieA 

whereas the matrix elements of the hopping operator are 

K n , n ' = -r) ^2 E <Wi'ei iCT ei:/ CT , ( 136 ) 
(ij)er<r=n 

where l ia = (0, . . . , 0, l ia , 0, . . . , 0) and © means mod 2 addition. By comparing (136) 
with (2) and observing that E(ij) 6 rE<7=t^n,n'®ii,ei^ = M> the unit value bein § 
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obtained only if n' E A(n), where 

A(n) = {n © \ ic © l ja : eT,a =|| with n ir7 © n jo = 1} (137) 
is the set of the configurations connected to n by the hopping of one particle, we have 

/ 1, n'eA(n), 

Vn, n . = V, = | Q) otherwige _ (138) 

In the following we will consider two dimensional square lattices having L x x L y sites 
and containing N p particles per spin. Periodic boundary conditions will be imposed. At 
the densities N p /(L x L y ) < 1/2 took into consideration, the set of all possible different 
values of the potential variables (7) during an infinitely long random walk corresponds 
to the set of all possible different values of V n , the matrix elements (135), over the 
configuration space. It is straightforward to see that Y = {0,7, 27, . . . , N p j}. The set 
of all possible different values of the hopping variables (17) during an infinitely long 
random walk, corresponds to the set of all possible different values of T n = rj\A(n) \ over 
the configuration space, with A(n) given by (137). The determination of 2? depends 
both on the number of particles and on the lattice size. For instance, in a lattice 2 x 3 we 
have & = {12r], I677, 20t?} with N p = 3 and 2? = {877, IO77, 12r], Urj, 16r]} with N p = 2. 
For hard-core bosons the set of all possible phase variables is ££ = {1}. 

Once we have determined the sets Y, 2? and , the asymptotic rescaled cumulants 
£( fc ) associated to considered Hubbard system can be measured as explained in [11]. 
Note that these cumulants are unaffected by a change of the parameters 77 and 7 of the 
Hamiltonian (134), whereas the label sets Y, 2? and ££ keep their cardinalities under 
the same change. These data are input into equations (127-128) to find the ground-state 
energy Eq. The results obtained by using cumulants up to order k max = 1,2,3,4 are 
shown in figure 2 as a function of the ratio 7/77 in the case of a lattice with 2x3 sites and 
N p — 2,3 particles per spin. In the same figure we also depict the values of E determined 
by exact numerical diagonalization. The curves obtained for /c max = 1 coincide with 
the uncorrelated multinomial prediction and, as already noted, their behavior is only 
qualitatively correct. Already at /c max = 2, the quantitative agreement with the exact 
ground-state energy becomes impressive at least for values of 7/77 not too large. By 
further increasing /c max , the quantitative agreement gradually improves in the whole 
range of 7/77 which, note the horizontal log scale, goes from the limit ||V"|| 1° 
the opposite one ||V"|| ^> In the case with N p — 2 particles per spin, the curve 

£0(7) obtained with /c max = 4 is indistinguishable from the reported exact values. In the 
case with N p — 3 particles per spin, the nonlinear equations (127-128) do not admit a 
solution with /c max = 3,4 when 7/77 is larger than a threshold value, namely 7/77 ~ 31.6 
for k max = 3 and 7/77 ~ 2.5 for A; max = 4. More precisely, above this threshold equations 
(127-128) admit only solutions with real negative or complex frequencies v. As discussed 
at the beginning of section 3, these are meaningless results which make the perturbative 
scheme invalid at the order considered. 

In figure 3 we show the results obtained with a system of larger size, namely a 
lattice 4x4 with N p = 2,4,5,8 particles. In this case, the number of configurations is 
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y/r\ 



Figure 2. Ground-state energy per particle for the 2x3 FNU hard-core boson 
Hubbard model with periodic boundary conditions versus the interaction strength 
7 with N p = 2 and N p = 3 particles per spin. We compare the results from 
exact numerical diagonalization ( x ) with those from present multinomial perturbative 
scheme by using cumulants up to order fc max = 1,2,3,4 (dotted, dashed, dot-dashed, 
solid lines, respectively). 



so large, namely 

^ Uv'J ' (139) 

that an exact numerical diagonalization of the Hamiltonian (134) is unfeasible. 
Therefore, we have compared the values of Eq predicted by the present multinomial 
perturbative scheme with those obtained by a Monte Carlo simulation [7]. The 
conclusions that we reach from the analysis of figure 3 are similar to those we noticed 
after figure 2. However, for this larger system we see that the parity of the order /c max 
may influence the quality of the approximation, a fact which is not surprising. For 
instance, in the case with N p = 5 it is evident that the curve £"0(7) obtained with 
fc max = 3 is not better than that obtained with /c max = 2. However, as evidenced in the 
enlargement shown in figure 4, the results obtained with /c max = 4 are more accurate 
than those with A; max = 2, at least in the 7/77 range where equations (127-128) admit 
the fc max = 4 solution. 

6. Conclusions 

In the framework of the probabilistic approach previously developed by us to study 
the ground-state properties of a many-body quantum system, we have introduced 
a multinomial perturbative scheme which has the following characteristics. At any 
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Figure 3. As in figure 1 in the case of a 4 x 4 FNU hard-core boson Hubbard model 
with periodic boundary conditions and N p — 2,4,5,8 particles per spin. The data 
indicated by o have been obtained by Monte Carlo simulations [7] (the associated 
statistical errors increase by increasing 7 and are of the order of the symbol size at 
7 ~ 1000?/). 




y/Ti 

Figure 4. Enlargement of figure 3, case with N p = 5. The perturbative scheme with 
^max — 4 admits a solution only for 7 < 3.2/7 but in this range provides results more 
accurate than those obtained with fc max = 1, 2, 3. 

perturbative order, the probability distribution of the potential, hopping and phase 
multiplicities, whose knowledge would allow for an exact solution of the problem, 
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is approximated by a multinomial-like distribution with infinitely many statistical 
moments. By increasing the perturbative order, an increasing number of cumulants 
of the distribution is made to coincide with the corresponding exact cumulants of the 
system. 

We have tested the proposed perturbative scheme in the case of two-dimensional 
lattice Hubbard models with spin l A hard-core bosons. Already at second perturbative 
order, we find a ground-state energy in good quantitative agreement with the exact one 
for any value of the ratio 7/77, 7 and r\ being the strengths of the interaction and hopping 
terms of the Hamiltonian of the system. The agreement improves at higher perturbative 
orders. At orders > 3, however, the scheme may not be always consistent, i.e. in some 
systems a solution for the ground-state energy is found only for values of 7/77 smaller 
than a threshold. As a matter of fact, we observe that in all our test cases at, or near, 
the Va particle filling, which is a case of remarkable physical interest, the solution of the 
perturbative method turns out to exist for all values of the coupling parameter up to 
the largest explored perturbative order, k max = 4, where it provides stunning results. 

The perturbative probability distribution at order k is built up from the knowledge 
of the first k connected statistical moments of the potential, hopping and phase 
multiplicities of the system. These cumulants are easily measured by Monte Carlo 
simulations and are independent of the parameters 7 and 77 of the Hamiltonian. Once 
the probability distribution is determined at the chosen approximation, the ground-state 
energy E , or any other correlation function in the ground state of the system, can be 
found by solving numerically a small system of nonlinear equations. This last job has 
a computational cost negligible with respect to the determination of the cumulants, 
which, in turn, has a cost equivalent to a direct Monte Carlo evaluation of E . Thus, 
the advantage of our approach in comparison to a direct Monte Carlo simulation is 
remarkable. Different Monte Carlo runs are needed to evaluate E for different values 
of 7 and/or 77, whereas in our approach we have to solve each time a small system of 
nonlinear equations and, una tantum, calculate the cumulants. 

The present perturbative probabilistic approach is particularly promising for 
systems affected by the so called sign problem, for which unbiased Monte Carlo 
simulations of E are impractical. In fact, the statistical evaluation of the cumulants of 
the potential, hopping and phase multiplicities is unaffected by the sign problem, which 
remains confined in the expression of the probability distribution and can be tackled by 
analytical techniques. We plan to discuss the case of fermions in a future paper. 

Appendix A. Solution of nonsymmetric algebraic Riccati equations 

In section 3.4 we have seen that the parameters p*- 2 * 1 , more precisely the associated 
reduced matrix p^ 2 \ are determined by the NARE (99). In general, the NAREs are 
defined as the quadratic matrix equations of the kind 



XCX - AX - XD + B = 0, 
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where we assume that the unknown X, as well as the coefficients A, B, C and D are 
quadratic matrices of finite size. In this section we illustrate two numerical methods 
developed to solve equation (A.l). The first one is an iterative method based on a fixed- 
point technique [22], whereas the second one is a direct method based on the Schiir 
decomposition [18]. 

Equations (A.l) play an important role in the study of the stochastic fluid models 
and have been extensively studied. In general, a NARE admits more than one solution. 
In most stochastic fluid models, the coefficients A, B, C and D form a super matrix 



H= t~> A ■ (A.2) 




with the property to be a so called M- matrix [18]. It can be proved that in this case 
the Schiir method provides the minimal non negative solution of the NARE, which is, 
there, the solution of physical interest. 

In our context, H is not a M- matrix and it is not clear which solution of the 
NARE (99) has to be considered. We propose to consider the unique solution given by 
the Schiir method. This solution coincides with that obtained by the iterative method 
in which X is chosen at the 0-th iteration as the solution of (A.l) with C = 0. Since in 
our case X represents the matrix of the perturbative parameters p^ 2 \ which we expect 
to be small, the above proposed solution seems the most natural one. 

Appendix A.l. Iterative method 

In [22] a class of fixed-point methods is considered to solve equation (A.l). These 
fixed-point iterations are based on a suitable splitting of the matrices A and D, that is 
A = Ai — A 2 and D = Di — D 2 , and have the form 

A ± X k+1 + X k+1 D ± = X k CX k + A 2 X k + X k D 2 + B, (A.3) 

with k — 0, 1, 2, . . . and X = 0. In particular, for A\ = A and D\ = D we have 

AX k+1 + X k+1 D = X k CX k + B, X = 0. (A.4) 

Note that finding X k+ i in terms of X k at k-th iteration implies to solve a Sylvester 
equation. This can be accomplished by vectorization, namely 

vec (AX k+1 ) + vec (X k+1 D) = vec (X k CX k + B) . (A.5) 

Using the properties of the vec operator, in particular 

vec (AXk+i) = (I <g) A) vec (X k+1 ) , (A.6) 

vec (X k+1 D) = (D T ® /) vec (X k+1 ) , (A.7) 

where <8> indicates the Kronecker product and T the transpose, equation (A.5) is 
rewritten as 

[(/ <g) A) + (D T ® /)] vec (X k+1 ) = vec (X k CX k + B) . (A.8) 
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This is a linear matrix equation which can be solved by standard methods, e.g. LU- 
factorization [19]. 

The convergence of the full class of iterative schemes (A. 3) to a solution X of 
(A.l) is ensured by a theorem [22]. In this class, the iterative scheme (A. 4) is the most 
expensive from a computational point of view, but, on the other hand, it has the highest 
(linear) convergence speed. 

Appendix A. 2. Schiir method 

In the following we discuss a different approach to solve equation (A.l) based on the 
ordered Schiir decomposition. This approach was conceived by Laub [17] for a symmetric 
algebraic Riccati equation and extended by Guo [23] to the study of NAREs. 
Let us rewrite the matrix H associated to the coefficients of (A.l) as 

Note that H is real in our case. We look for an orthogonal transformation 

(A10) 

which leaves H in a semi-ordered real Schiir form, 

U T HU = S = ( Sn f 12 V (A.ll) 



S- 



22 



in which S n and S22 contain only blocks, denoted s^-, i,j = 1,2,..., of size 1 or 2. The 
eigenvalues of the 2x2 diagonal blocks 8 a provide the complex conjugated eigenvalues 
of H whereas the lxl blocks are the real eigenvalues of H . The diagonal blocks are 
semi-ordered in the sense that if Su, Sjj and Sj~k have eigenvalues with positive, null and 
negative real parts, respectively, then i < j < k. It is possible to show that the matrix 
U n is invertible [24] and that 

X = U 2 iU n 1 (A. 12) 

solves (A.l). Note that the semi-ordered decomposition (A.ll) is unique and so is 
the solution (A. 12). We used the subroutines of LAPACK library [25] to numerically 
implement the Schiir method. 

Appendix B. Equations for the perturbative parameters: fourth order 

The perturbative parameters are determined by the system of equations 

((^a 2 ^ a MP (1 \p {2 \p {3 \p {4) ) = (B.l) 

with ai, ai2, CK3, 0:4 G Jif. By using (73) and taking the derivative of (80) with respect 
to J ai , the above system can be cast in the form 

V V V V fA (2 ' 0) £ {2) S (2) S (2) +£ (2 ' 0) A {2) S (2) S {2) 

a£,jr peje jeje sej? 
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, E (2,0) E (2) a(2) v,(2) v(2,0)v( 2 ) Yj(2) a(2) \ (4) _ a (DO) 

where E^ 2 ) is the asymptotic rescaled cumulant of order 2 and the matrices S( 2 '°), 
A^ 2 ^ and A^ 2,0 ) are defined by (90), (104) and (105), respectively. The tensor A has 
components a±, a 2 , a 3 , a 4 G Jf given by 

A = V(4) _ W4,0) 

^ J oiQ2"3 Q! 4 /J o 1 a2a3a4 /J oia2a3a4 



where 



-(2,0) / M (3) v.(3) ^(3) „(3) 



- V V V r( 2 <°) F (3) r {2) r( 3 ) 

+ S^FW jE « E« Eg 4 , (B.3) 

ae^r /Je^r 7£jT <5e 

(4,0)_ y (3) A / ( 3 ) v (3) v (l) (2W2) (2) (2)\ , R ,x 

(2) (2) (2) (2) 

7 6JT ^7 



r(3 ) _ ^(^P) 



(3) (3) (3) (2) (2) (2) (2) (2) (2) 

F (3) _ Pa/3-/ Pafij Pa(3j _ PgpP°n _ PafiPfr _ P<*~lP~,p 

**fh~ y(l) + y(l) + y(l) ' (1)2 (1) 2 (1) 2 

^7 La ^7 

/ (3) (2) (3) (2) (3) (2) (2) (2) (2) \ 

El Pafi&PS"/ . PaS-yPsp . Psp-yPaS n Pa5 Ps/3 Ps-y \ /T3 „\ 

v^r t " 2 ^r ) • ( ' 

(3) (2) . (3) (2) . (3) (2) (3) (2) . (3) (2) . (3) (2) 

p(4) _ PaPjPaS PayPal ^ Pg^sPaP PaPjP/35 ~^ PafiSPpf ~^ PpjSPqfi 

a ^ S ~ v (l) 2 + v (l) 2 

(3) (2) . (3) (2) . (3) (2) (3) (2) . (3) (2) . (3) (2) 

Pg^PjS ^ Pg-f&P-tP PpjSPai PafjjPs-y ~T PgSjPsp Psp-yPaS 

+ ^(D 2 + ^(D 2 

^7 ^6 

(2) (2) (2) (2) (2) (2) (2) (2) (2) (2) (2) (2)' 

P apPaiPuS Paf3Pf3-/Pl35 PaiP-yfjP^S PaSPsp P5j 

y(l) 3 ^(l) 3 Ul) 3 Ul) 3 

Zjq, Zj^j Zj^y Zj£ 
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■ (3) (3) (3) (3) (3) (3) 

PafiePe-fS + Pa£-yP £ /3S + P £ /3~fPea5 



(3) (2) (2) . (3) (2) (2) . (3) (2) (2) 
PaPePe-rPeS + P^PJpJ + P £ frP<*eP e 5 

(3) (2) (2) . (3) (2) (2) . (3) (2) (2) ' 
PaeSPepPzi + Pjfpeclph' + P^pUp^ 



+ 6 



+ 



(2) (2) (2) (2) 

s^ 3 



(B.8) 



To find p^\ we first determine the reduced tensor p^ which is the solution of the linear 

system 

V V V V (\(my( 2 ) V(2) f(2) ,y(2,0)T(2) y(2) y( 
Z^ Z^ Z^ Z^ l/Va ^l3a2^-ya 3 ^Sa 4 + ^aia A /3<*2 ^s 2 ^ 

a£ J# /36 J# 76 J# <5e 



^(2) 
J <5a 4 



, ^(2,0)f (2) 7(2) y(2) ^(2,0)y(2) y(2) T( 



(2,U)y^j yj(2j a (2) \ 44) _ a / r q n 



with a 1 ,a 2 ,a 3 ,a 4 e The matrices £ (2) , £ (2 - 0) , A( 2 ) and A (2,0) are defined by (94), 
(96), (108) and (109), respectively. The complete fourth order perturbative parameter 
is then recovered using the sum rules (88) for k — 4. 
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